#------------------------------------------------------------------------------
rm(list = ls())
library(LalRUtils)
libreq(
  data.table, zoo, tictoc, fixest, PanelMatch, patchwork,
  rio, magrittr, janitor, did, panelView, ggplot2, RPushbullet, ggiplot, 
  tidyverse, data.table, zoo, tictoc, fst, fixest, PanelMatch, patchwork,
  rio, magrittr, janitor, did, panelView, ggiplot, tictoc, binsreg, interflex
)

set.seed(42)
theme_set(lal_plot_theme())

notif = \(x) pbPost("note", x)

#------------------------------------------------------------------------------




#------------------------------------------------------------------------------
# Define Paths
#------------------------------------------------------------------------------

# R studio
setwd( dirname(rstudioapi::getActiveDocumentContext()$path) )
# R default : unccoment if you use default R
# setwd(getSrcDirectory(function(){})[1])
#------------------------------------------------------------------------------



#------------------------------------------------------------------------------
# Load Data
#------------------------------------------------------------------------------

vcf <- fread("vcf_data_complete.csv", sep = ",")
setnames( vcf, "d", "D")


#------------------------------------------------------------------------------



#------------------------------------------------------------------------------
# Figure 5: Dynamic Treatment Effects of PESA Adoption on Forest Index
#------------------------------------------------------------------------------

if (exists("vcf_data")) {
  vcf = vcf_data[year >= 1995]
  rm(vcf_data)
}

# %% GFC
vcf[, time := year - first_pesa_exposure]
vcf[, pref_bin := ntile(cover_1990, 10)]
vcf_evsamp = vcf[time %between% c(-6, 6)]
# %%
estudy_plot = function(cutoff, title){
  es1 = feols(forest_index ~ i(time, sch,  ref=-1) | cellid + styear,
              cluster = ~ blk, vcf_evsamp[pref_bin >= cutoff])
  f = ggiplot(es1) + ylim(c(-3, 3)) + ggtitle(title) + lal_plot_theme()
  return(f)
}
# %%
ff = (estudy_plot(1, "Full sample") + labs(y = "", x = "") + 
        theme(axis.text = element_text(size = 20), 
              plot.title = element_text(size = 22))|
        estudy_plot(6,  "5th decile and up") + labs(y = "", x = "") + 
  theme(axis.text = element_text(size = 20), 
        plot.title = element_text(size = 22)) ) /
  (estudy_plot(7,  "6th decile and up") + labs(y = "", x = "") + 
     theme(axis.text = element_text(size = 20), 
           plot.title = element_text(size = 22)) |
     estudy_plot(8,  "7th decile and up") + labs(y = "", x = "")  + 
  theme(axis.text = element_text(size = 20), 
        plot.title = element_text(size = 22)) ) /
  (estudy_plot(9,  "8th decile and up") + labs(y = "", , x = "Time") + 
     theme(axis.text = element_text(size = 20),
           plot.title = element_text(size = 22), 
           text = element_text(size = 20)) |
     estudy_plot(10, "9th decile and up") + labs(y = "", x = "Time") + 
  theme(axis.text = element_text(size = 20), 
        plot.title = element_text(size = 22), 
        text = element_text(size = 20)) )
ggsave("main_figure5.pdf", ff, width = 12, height = 15,
       device = cairo_pdf)

#------------------------------------------------------------------------------
